
# 
# This script reproduces and applies code from the replication files of Yimeng Li (2019)
# "Relaxing the No Liars Assumption in List Experiment Analyses"
# DOI: https://doi.org/10.1017/pan.2019.7

remove(list=ls())

# Load libraries:
library(bindata) # required by sim_bounds.R
library(dplyr)
library(forcats)
library(haven)
library(ggplot2)
library(list)
library(here)
library(nleqslv) # required by list_relaxed_liars.R
library(purrr)
library(tidyr)
library(patchwork)

source("Li_2019_NoLiars/list_relaxed_liars.R")
source("Li_2019_NoLiars/list_types_plot_DE.R")
source("Li_2019_NoLiars/list_create_df.R")


list_dataDE <- read_dta("data/listDE.dta")
list_dataDE <-as.data.frame(list_dataDE)
list_dataDE<- list_dataDE  %>% 
  mutate(treatfac=as.factor(treated),
         IDvar=as.factor(IDvar),
         round.fac=as.factor(round)) |> 
  mutate(value= value-1)

options(digits = 3)


dist_control_Example <- c(419, 748, 146, 31)
dist_treated_X_Example <- c(182, 439, 594, 113, 19)

DiM_X_Example <- ictreg(y ~ 1, data = create_list_df(dist_control_Example, dist_treated_X_Example, J = 3), treat = "treat", J = 3, method = "lm")

unname(DiM_X_Example$par.treat)
unname(c(DiM_X_Example$par.treat+qnorm(0.025)*DiM_X_Example$se.treat, DiM_X_Example$par.treat+qnorm(0.975)*DiM_X_Example$se.treat))

bounds_X_Example <- list_relaxed_liars(dist_control_Example, dist_treated_X_Example, J = 3, affirmative = TRUE)

bounds_X_Example$bounds
bounds_X_Example$CI

FigureDE<- list_types_DE_plot(bounds_X_Example) 

F1<- FigureDE$plot_two_types
F2<- FigureDE$plot_three_types

F1+F2
ggsave("FigureA8.png",
       path = here("figures"),
       dpi=600)
